home *** CD-ROM | disk | FTP | other *** search
/ CD-ROM Today - The Disc! 5 / CD-ROM Today - The Disc (Issue 5)(November 1994).ISO / mac / Mac shareware / Education / RLaB / toolbox / pascal.r < prev    next >
Encoding:
Text File  |  1994-09-21  |  1.2 KB  |  52 lines  |  [TEXT/ttxt]
  1. //
  2. //PASCAL  PASCAL(N) is the Pascal matrix of order N: a symmetric positive
  3. //        definite matrix with integer entries, made up from Pascal's
  4. //        triangle.  Its inverse has integer entries.
  5. //        [Conjecture by NJH: the Pascal matrix is totally positive.]
  6. //        PASCAL(N,1) is the lower triangular Cholesky factor (up to signs
  7. //        of columns) of the Pascal matrix.   It is involutary (is its own
  8. //        inverse).
  9. //        PASCAL(N,2) is a transposed and permuted version of PASCAL(N,1)
  10. //        which is a cube root of the identity.
  11.  
  12. //    J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  13. //    Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
  14.  
  15. pascal = function(n, k)
  16. {
  17.   local(P, i, j);
  18.  
  19.   if (!exist (k)) { k = 0; }
  20.  
  21.   P = diag( (-1).^[0:n-1] );
  22.   P[;1] = ones(n,1);
  23.  
  24.   // Generate the Pascal Cholesky factor (up to signs).
  25.   for( j in 2:n-1 )
  26.   {
  27.     for( i in j+1:n )
  28.     {
  29.       P[i;j] = P[i-1;j] - P[i-1;j-1];
  30.     }
  31.   }
  32.  
  33.   if( k == 0 ) 
  34.   {
  35.     P = P*P';
  36.   else if( k == 2 ) {
  37.     P = P';
  38.     P = P[ n:1:-1 ;];
  39.     for( i in 1:n-1 )
  40.     {
  41.       P[i;] = -P[i;];
  42.       P[;i] = -P[;i];
  43.     }
  44.     if( n/2 == round(n/2) ) 
  45.     {
  46.       P = -P;
  47.     }
  48.   }}
  49.  
  50.   return P;
  51. };
  52.